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Abstract 

This  paper  uses  KAM  torus  theory  and  Simplified  General  Perturbations  4  (SGP4) 
orbit  prediction  techniques  compiled  by  Dr.  William  Wiesel  and  compares  it  to  Ana¬ 
lytical  Graphics®  Incorporated  (AGI)  Satellite  Toolkit®  (STK)  orbit  data.  The  goal 
of  this  paper  is  to  verify  KAM  torus  theory  can  be  used  to  describe  and  propagate  an 
Earth  satellite  orbit  with  similar  accuracy  to  existing  general  perturbation  techniques. 
Using  SGP4  code  including  only  truncated  geopotential  effects,  KAM  torus  generating 
code,  and  other  utilities  were  used  to  describe  a  particular  satellite  orbit  as  a  torus  and 
then  propagate  the  satellite  using  traditional  and  KAM  torus  techniques.  These  results 
were  compared  with  similar  data  generated  from  initial  conditions  in  STK.  Comparisons 
show  orbit  prediction  for  this  particular  satellite  can  be  made  with  low  kilometer  level 
accuracy.  It  is  claimed  with  increased  mathematical  precision  and  orbital  model  detail, 
KAM  torus  theory  applied  to  orbit  prediction  techniques  can  produce  more  accurate 
results  than  currently  achievable. 
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VERIFICATION  OF  KAM  THEORY 
ON  EARTH  ORBITING  SATELLITES 

I.  Introduction 

The  Air  Force’s  Joint  Space  Operations  Center  is  tasked  with  constantly  providing 
United  States  Strategic  Command  a  space  situational  awareness  picture.  All  objects 
orbiting  the  Earth  ten  centimeters  across  or  larger  are  detected,  tracked,  identified,  and 
cataloged.  The  sum  total  of  all  these  objects  is  a  computerized  library  of  over  17,000 
man-made  objects  called  the  Space  Catalog.  JSpOC  personnel  rely  on  only  twenty-nine 
space  surveillance  sensors  worldwide  to  accurately  maintain  position  data  for  all  objects 
in  the  Space  Catalog  [3] .  The  sheer  number  of  items  requiring  tracking  with  so  few  assets 
shows  how  daunting  a  task  this  can  be. 

1.1  Motivation 

Current  general  perturbation  orbit  prediction  techniques  produce  acceptable  results 
for  low  priority  objects.  However,  there  is  plenty  of  room  for  improvement.  Higher 
priority  missions  such  as  manned  space  walks,  collision  avoidance,  and  atmospheric  re¬ 
entry  already  require  numerical  integration  for  orbit  prediction  [4],  In  either  case,  as  time 
increases  so  does  location  uncertainty.  The  location  uncertainty  is  an  ellipsoid  in  which 
the  object  most  likely  resides.  Location  determination  techniques  produce  an  estimate  of 
an  objects  location.  Assumptions  or  errors  residing  in  the  technique  or  equipment  used 
correspond  to  location  errors.  In  Figure  1.1,  the  location  uncertainty  (blue  ellipsoid)  of  a 
reference  satellite  is  shown  growing  as  time  increases.  The  object  must  be  detected  and 
tracked  again  before  the  ellipsoid  grows  too  large  and  the  satellite  is  lost.  If  a  satellite 
was  lost,  meaning  its  uncertainty  ellipsoid  had  grown  too  large,  resources  and  time  would 
be  required  to  meticulously  scan  the  entire  ellipsoid  to  reacquire  the  satellite. 

Spacetrack,  a  phased  array  radar,  is  used  to  track  objects  in  space.  A  quick 
search  for  “Spacetrack  problems”  on  the  Internet  will  show  multiple  articles  detailing 
new  scheduling  methods  and  algorithms  to  track  as  many  high  priority  objects  as  possi- 
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ble  within  a  given  time  period.  This  research  is  being  completed  since  tracking  all  objects 
with  the  precision  desired  is  becoming  impossible  [5].  These  articles  show  current  space 
object  tracking  is  already  a  problem  or  soon  will  be. 


(e)  (f) 

Figure  1.1:  Reference  satellite  error  ellipsoid  increasing  with  time. 


On  February  10,  2009,  a  Russian  satellite,  Cosmos  2251,  collided  with  a  U.S. 
communications  satellite,  Iridium  33,  in  orbit.  This  event  produced  thousands  of  pieces 
of  space  junk  which  required  tracking.  It  also  increased  the  chances  of  one  of  the  other 
36  Iridium  satellites  in  similar  orbits  being  hit  by  50%  [6].  A  simulated  picture  moments 
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after  the  collision  can  be  seen  in  Figure  1.2.  Events  such  as  this  set  to  quickly  overload 
current  tracking  resources. 


Predicted  debris  trajectories 

Statistical  model  (Gaussian  distribution) 

Aagi 

Debris  not  to  scale 

Initial  Iridium  33 
Trajectory 

1  km/s  delta-V 

\  ^ 

Az  -90  deg,  El  +/-  38  deg 

Initial  Cosmos  2251 

JGM2  (12x12)  gravity  field 

T  rajectory  \ 

Figure  1.2:  A  model  showing  approximate  aftermath  of  the  Cosmos  /  Iridium  collision. 
“Image  courtesy  of  Analytical  Graphics,  Inc.  (www.agi.com).” 


1.2  Problem  Statement 

One  option  for  the  overloaded  tracking  situation  is  to  build  more  tracking  sensors 
and  employ  more  people  to  operate  and  maintain  them.  This  option  would  cost  large 
amounts  of  money  on  both  the  equipment  and  manning  side.  Additional  computer 
resources  would  also  need  to  be  procured.  Another  option  would  be  to  work  smarter 
and  find  a  way  to  accurately  predict  these  orbiting  object’s  position  further  out  into 
the  future,  lessening  the  frequency  with  which  their  tracks  need  to  be  updated.  Using 
a  theory  conceived  almost  60  years  ago,  Dr.  William  Wiesel  at  the  Air  Force  Institute 
of  Technology  may  have  found  such  a  method.  It  utilizes  Kolmogorov-Arnol’d-Moser 
(KAM)  theory  which  is  named  after  the  discoverer,  A.  N.  Kolmogorov,  and  subsequent 
verifiers,  V.  I.  Arnol’d  and  J.  Moser.  Initial  work  by  Dr.  Wiesel  shows  Earth  satellite 
orbits  can  be  described  using  KAM  theory  [7].  The  current  research  attempts  to  validate 
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the  theory  against  simulated  ‘real’  satellite  data  produced  using  Analytical  Graphics®, 
Incorporated  (AGI)  Satellite  Toolkit®  (STK). 

1.3  Method  of  Investigation 

The  object  of  this  paper  is  to  test  satellite  prediction  using  KAM  theory  against 
current  methods.  The  general  perturbations  method  called  Simplified  General  Pertur¬ 
bations  4  (SGP4)  will  be  used  to  predict  satellite  location  as  the  current  method.  Given 
the  algorithm  necessary  for  SGP4  to  produce  accurate  data  for  intervals  greater  than  a 
few  days,  SGP4  and  the  KAM  theory  method  cannot  be  directly  compared.  To  solve  this 
problem  STK  will  be  used  as  an  intermediary.  STK  will  be  compared  to  SGP4  to  ensure 
similar  or  better  results  are  achieved  with  the  more  accurate  STK  numerical  integrator. 
STK  orbit  data  will  then  be  compared  to  the  KAM  theory  method.  If  similar  results  are 
achieved,  then  it  will  have  been  shown  that  KAM  theory  orbit  prediction  is  at  least  as 
accurate  as  the  current  SGP4.  As  an  initial  step  it  is  important  to  verify  KAM  theory 
can  be  applied  to  Earth  orbiting  satellites.  Once  this  is  complete  more  detail  can  be 
added  to  the  KAM  theory  method  to  produced  better  results  than  currently  achievable 
in  less  time. 

Actual  Earth  satellites  will  be  used,  but  their  location  data  will  be  produced  from 
AGI’s  STK  software.  Previous  investigations  by  Little  [8]  used  actual  satellite  location 
data  provided  from  a  satellite  team,  but  problems  were  encountered  when  satellite  ma¬ 
neuvers  were  made.  Using  the  STK  software,  accurate  location  data  can  be  used  while 
removing  any  outside  influence  on  the  problem,  allowing  the  two  methods  to  be  more 
accurately  compared. 

1.4  Results 

Comparison  between  KAM  torus  orbit  prediction  techniques  and  STK  data  as  a 
truth  model  resulted  in  position  differences  under  3.5  kilometers  out  to  120  days.  This 
research  is  another  source  that  proves  Earth  orbiting  satellites  can  be  modeled  as  a  torus 
using  KAM  theory  and  propagated  accurately.  The  next  step  would  be  to  incorporate 
more  detail  into  the  KAM  torus  model  and  compare  to  similar  detailed  STK  data.  Upon 
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successful  comparison,  the  detailed  KAM  torus  model  would  be  compared  with  an  actual 
satellite  and  its  positional  data.  Given  the  accuracy  level  desired,  quadruple  precision  is 
most  likely  required. 
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II.  Background 


2.1  Hamiltonian  Dynamics 

From  the  publication  of  Sir  Isaac  Newton’s  Philosophiae  Naturalis  Principia  Mathe- 
matica  in  1687  np  until  the  late  eighteenth  century,  Newtonian  Mechanics  was  mechan¬ 
ics  [9].  With  F  =  ma,  forces  on  single  particles,  or  systems  of  particles,  were  used 
to  solve  problems  of  the  day.  However,  Newtonian  mechanics  had  its  drawbacks  and 
some  problems  could  not  be  solved  easily.  In  1788,  Joseph-Louis  Lagrange  published 
Mecanique  Analytique  and  analytical  mechanics  was  born.  In  1834,  William  Hamilton 
released  two  papers  on  general  methods  in  dynamics  in  which  he  applied  his  principle 
to  the  Lagrangian  [9].  This  was  manipulated  into  the  Hamiltonian  and  Hamiltonian’s 
equations  of  motion,  which  described  the  system  [10]. 

It  should  be  noted  that  in  this  paper,  systems  are  assumed  to  be  holonomic,  con¬ 
servative  systems.  The  material  in  this  paper  does  apply  to  other  systems  as  well,  but 
this  constraint  helps  lessen  the  complexity.  A  few  topics  must  be  covered  before  jumping 
right  into  the  Hamiltonian,  the  first  being  generalized  coordinates. 

2.1.1  Generalized  Coordinates.  Most  people  are  familiar  with  the  standard 
Cartesian  coordinate  system  of  x,  y,  and  z.  However,  there  are  many  other  types  of 
coordinate  systems:  polar,  spherical,  and  cylindrical  to  mention  a  few.  These  coordinate 
systems  help  describe  the  position  of  an  object  in  a  particular  frame  of  reference  or  space. 
Generalized  coordinates  can  be  thought  of  as  a  way  of  describing  a  system  without  being 
restricted  to  physical  positional  coordinates.  This  idea  expanded  how  a  system  could  be 
described  and  allowed  problems  to  be  solved  in  many  new  ways. 

The  coordinates  of  a  system  can  indeed  be  the  x,  y,  and  z  position  of  an  object 
in  space,  but  they  could  also  be  the  frequencies  describing  its  periodic  motion,  the 
coefficients  of  a  Fourier  series  expansion,  a  set  of  angles,  or  almost  anything  else.  The 
concept  of  generalized  coordinates  allows  the  coordinates  to  be  almost  anything.  Most 
problems  can  be  solved  in  multiple  coordinate  systems  and  using  a  different  system  may 
make  a  problem  much  easier  to  solve.  Changing  from  physical  to  non-physical  coordiantes 
could  give  new  insight  to  an  existing  problem  or  make  its  solution  easier.  Generalized 
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coordinates,  q,  are  most  commonly  referred  to  as, 


9  =  (9i,92,93,-..,9„), 


(2.1) 


and  for  the  problem  presented  here  are, 


(  \ 

x 


q  =  S  y 


z 

\  / 


(2.2) 


The  number  of  coordinates  helps  determine  the  dimensionality  of  a  problem. 


2.1.2  Dimensionality.  Another  concept  that  must  be  grasped  is  dimensionality. 
We  all  live,  work,  and  are  most  familiar  in  three-dimensional  space,  ffowever,  mathe¬ 
matical  problems  are  not  restricted  to  three  dimensions.  A  simple  pendulum  could  have 
three  coordinates  while  a  more  complicated  problem  could  have  six  coordinates  or  more. 
Although  a  six  dimensional  object  cannot  be  visualized  in  three  dimensional  space,  the 
mathematics  behind  the  problem  work  just  fine.  Later  in  this  paper  when  six  dimensional 
objects  are  discussed,  think  of  a  three  dimensional  object  for  visualization  purposes,  but 
realize  the  actual  dimensionality  is  greater. 


2.1.3  The  Lagrangian.  As  mentioned  earlier,  Lagrange  initiated  what  is  now 
know  as  analytical  mechanics,  ffowever,  he  relied  on  Jean  le  Rond  d’Alembert’s  1743 
work,  “Traite  de  dynamique” ,  where  d’Alembert  revealed  his  principle, 

N 

y:  (Fk  -  pk)  .  5rk  =  0.  (2.3) 

k= 1 

With  this  principle,  a  dynamic  problem  is  treated  as  a  static  one.  The  total  virtual  work 
done  on  a  system  by  infinitesimal  virtual  displacements  is  equal  to  zero  [10]. 

Using  d’Alembert’s  principle  and  work  by  Euler,  Lagrange  began  thinking  of  the 
system  as  a  whole  and  the  total  energy  of  the  system,  instead  of  individual  particles. 
Energy  conies  in  two  forms,  potential  and  kinetic.  Kinetic  energy  (T)  is  the  energy  of  a 
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system  based  on  its  motion  and  is  given  by, 


X 

T  =  -mv2.  (2.4) 

While  potential  energy  ( V )  is  the  energy  stored  within  the  system.  For  systems  near  the 
surface  of  the  Earth,  the  potential  energy  will  be  based  on  position  and  can  be  given  by, 

V  =  mgh.  (2.5) 


However,  in  dealing  with  celestial  bodies  and  their  gravitational  fields  a  different  approach 
must  be  used.  For  the  moment  the  above  representation  is  sufficient,  but  a  more  accurate 
model  will  be  discussed  in  Section  2.2.1. 

A  cornerstone  to  analytical  dynamics  is  a  mathematical  expression  called  the  La- 
grangian  (<£),  which  is  the  difference  between  the  kinetic  and  potential  energy  of  a 
system, 

£  —  T  —  V.  (2.6) 


The  Lagrangian  is  called  a  descriptive  function  for  a  system,  because  it  contains  all  the 
information  necessary  to  produce  the  equations  of  motion  for  the  system.  However,  the 
Lagrangian  produces  n  second-order  differential  equations,  where  n  is  the  number  of 
coordinates  for  the  system, 


d  fd£ 
dt  \dqn 


d£ 

dqr 


=  0. 


(2.7) 


While  second-order  differential  equations  can  be  solved,  first-order  differential  equations 
would  be  easier  to  solve. 


The  time  derivatives  of  generalized  coordinates  are  generalized  velocities,  q.  The 
derivative  of  the  Lagrangian  with  respect  to  the  generalized  velocities  are  the  generalized 
momenta, 


Pn  = 


d£ 

d4n 


(2.8) 


Generalized  coordinates  and  generalized  momenta  come  in  pairs,  thus  they  are  normally 
referred  to  as  conjugates  of  each  other.  If  a  generalized  coordinate  or  momenta  does  not 


appear  in  the  Lagrangian,  then  its  conjugate  is  cyclic,  or  ignorable,  and  is  constant.  If  a 
generalized  coordinate  or  momenta  is  constant,  a  conservation  law  is  within  the  problem 
and  much  can  usually  be  determined  about  the  system  without  ever  having  to  solve  the 
problem  [10]. 

2.1.4  Hamilton’s  Principle.  William  Rowan  Hamilton  began  applying  varia¬ 
tional  calculus  to  previous  work  in  the  emerging  field  of  analytical  mechanics.  He  looked 
at  the  entire  dynamical  system  between  two  points  in  time,  t±  and  t2  [10].  With  varia- 


Figure  2.1:  Hamilton’s  Principle 


tional  calculus,  the  path  between  these  two  points  can  vary  as  long  as  the  endpoints  stay 
fixed.  Integrating  over  this  distance  produces  a  scalar  integral  which  is  independent  of 
the  coordinates  of  the  Lagrangian  [10], 


(2.9) 


£dt  =  0. 


Hamilton’s  principle  is  essentially  an  integrated  form  of  d’Alembert’s  Principle.  Instead 


of  using  d’Alembert’s  Principle  to  derive  Lagrange’s  equations  of  motion,  Hamilton’s 
Principle  offered  an  easier  way  which  was  not  dependent  on  physical  coordinates  [10]. 

2.1.5  Hamilton’s  Equation.  While  Lagrange’s  equations  of  motion  are  a  set 


of  n  second-order  differential  equations,  Hamilton’s  equations  of  motion  are  a  set  of  2 n 


first-order  differential  equations  [10].  The  Hamiltonian,  PL,  can  initially  be  defined  as, 


(2.10) 
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To  get  first-order  differential  equations  and  the  Hamiltonian  in  the  correct  form,  the 
conjugate  momenta,  Equation  2.8,  must  be  substituted  into  Equation  2.10  [10], 


hi  =  ^PkQk  -  £■  (2.11) 

fc=i 


The  Hamiltonian  equations  of  motion  are  then  the  partial  derivatives  of  the  Hamiltonian 
with  respect  to  the  generalized  momenta  and  negative  generalized  coordinates  [10], 


qk  = 


m 

dpi’ 


(2.12) 


*  =  -ai-  <2-13) 

Finally  a  set  of  2 n  first-order  differential  equations,  also  known  as  Hamilton  canonical 
equations,  are  ready  to  be  solved.  The  combination  of  variables  qk  and  pk  that  solve  the 
above  equations  can  be  plotted,  and  the  trajectory  of  these  points  through  time  shows 
the  motion  of  the  system.  These  points  can  then  be  plotted  in  phase  space  [10]. 


2.1.6  Phase  Space.  For  mathematical  problems  in  grade  school,  Euclidean 
space  was  most  likely  used  because  its  properties  were  well  suited  for  geometrically 
showing  solutions  to  the  mathematics  being  worked  on  at  the  time.  When  dealing  with 
Hamiltonian  mechanics,  phase  space  is  used.  The  spot  in  phase  space  where  the  general¬ 
ized  coordinates  and  conjugate  momenta  produce  a  solution  to  the  Hamilton  canonical 
equations  a  point  is  produced.  The  combination  of  all  these  points  represents  the  mo¬ 
tion  of  the  system  in  phase  space.  While  Euclidean  space  has  three  dimensions  when 
n  equals  three,  phase  space  would  have  2 n  dimensions,  or  six.  Since  graphing  in  the 
six-dimensional  space  is  impossible,  phase  space  plots  usually  include  only  two  or  three 
of  the  variables  at  a  time  [10]. 


2.2  Earth  Orbiting  Satellite 

2.2.1  The  Geopotential.  The  geopotential  for  a  given  celestial  object  is  its 
gravitational  held.  This  gravitational  held  exerts  a  force  on  objects  based  on  their 
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distance  from  the  celestial  body  and  the  mass  of  the  celestial  body.  An  example  of  the 
two  body  problem  is  a  small  object  orbiting  around  a  larger  object.  The  small  object,  a 


Figure  2.2:  The  Two  Body  Problem 


satellite  for  example,  is  a  distance  r  from  the  center  of  the  larger  object,  the  Earth  for 
example.  The  satellite  is  traveling  with  a  velocity,  v,  and  should  continue  in  a  straight 
line  off  the  page.  However,  the  mass  of  the  Earth  creates  a  gravitational  field  in  which 
an  attractive  force  is  exerted  upon  the  satellite.  This  attractive  force  and  the  velocity  of 
the  satellite  form  an  equilibrium  that  allows  the  satellite  to  stay  in  an  orbit  around  the 
Earth.  The  attractive  force  the  Earth  exerts  on  the  satellite  is  the  potential  energy  and 
for  the  two  body  problem  is  given  by, 


V  —  (2.14) 

r 

where  r  is  the  distance  from  the  satellite  to  the  center  of  the  Earth  and  /i  is  a  gravitational 
parameter  for  the  Earth,  which  is  based  on  the  mass  and  gravitational  constant  of  the 
Earth  [10]. 

The  two  body  gravity  model  provides  decent  estimates  for  simplified  problems. 
However,  more  complicated  problems  require  a  better  estimate  for  the  geopotential.  In 
this  research  the  geopotential  will  be  estimated  using  the  following  equation, 

oo  n 

V  =  --  S)(Cnm  cosmA  +  S^sinmA),  (2.15) 
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where  the  Earth’s  gravitational  parameter  is  again  /j  and  the  Earth’s  equatorial  radius  is 
R®.  The  Legendre  polynomials,  P™,  and  the  held  coefficients,  Cnm  and  Snrn)  are  used  to 
help  describe  the  geopotential  more  precisely  anywhere  around  the  Earth’s  sphere.  The 
radius  to  the  satellite  from  the  center  of  the  Earth  (r),  the  geocentric  latitude  (5),  and 
the  east  longitude  (A)  are  found  from  [7], 

r  =  \f  x2  +  y2  +  z2,  (2-16) 

sin5=^=^’  (2-17) 

V  x2  +  y2 

tanA=— .  (2-18) 

x 

The  geocentric  latitude  is  the  angle  between  the  equatorial  plane  and  a  line  from  the 
center  of  the  Earth  to  a  point  on  the  surface.  This  value  is  less  than  the  geodetic  latitude 
more  commonly  used  which  is  the  angle  between  the  equatorial  plane  and  a  line  normal 
to  the  Earth’s  surface. 


The  gravity  model  used  by  Dr.  Wiesel’s  SGP4  is  EGM96,  or  the  Earth  Gravity 
Model  completed  in  1996.  For  Equation  2.15,  terms  to  degree  (m=20)  and  order  (n=20) 
are  used.  Although  sufficient  for  application  in  this  research,  the  EGM96  model  does  go 
up  to  order  and  degree  360.  An  example  of  the  variation  in  geopotential  can  be  seen  in 
Figure  2.4  [2], 
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Figure  2.4:  GRACE  Gravity  Model,  Europe  and  Africa  [2] 


2.2.2  Hamiltonian.  For  an  Earth  satellite  Hamiltonian,  an  Earth-Centered, 
Earth-Fixed  (ECEF)  reference  frame  will  be  used.  Earth-Centered  means  the  center 
of  the  Earth  is  at  the  origin  of  the  reference  frame,  and  Earth-Fixed  means  the  Earth 
appears  stationary  with  respect  to  the  reference  frame.  In  actuality  the  Earth  is  still 
revolving,  however,  since  the  reference  frame  is  revolving  with  the  Earth  it  looks  fixed. 
The  Hamiltonian  can  be  written  as  [7], 

ft  =  \(p2x  +  P2y  +  P2z)  +  ve(ypx  -  xpy ) 

00  71  /  \  —n 

-  -  V!  (  jr  )  Pn  (sin  <5)  ( Cnm  cos  m\  +  Snm  sin  m\) .  (2.19) 

The  Earth’s  rotational  angular  velocity  is  given  by  0%.  The  entire  last  term  of  Equation 
2.19  is  the  Earth’s  geopotential  as  described  earlier.  The  generalized  coordinates  and 
momenta  are  then, 


'  ' 

X 

/ 

Px 

q  =  < 

y 

> , 

(2.20) 

p  =  < 

Py 

z 

Pz 

\ 
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As  was  the  Lagrangian,  the  Hamiltonian  is  also  a  descriptive  function  in  that  the 
equations  of  motion  for  the  system  can  be  obtained  from  it.  These  equations  of  motion 
are  called  Hamilton  canonical  equations  [10], 


/ 

—  CJ^Py  ~ 

dV 

dx 

/ 

Px 

+  u®y 

U®Px  ~ 

8V 

dy 

>  , 

(2.22) 

p  =  < 

Py 

—  (J&X 

dV 

Pz 

k  dz 

J 

< 

(2.23) 


A  complete  derivation  of  the  above  Hamiltonian  and  canonical  equations  can  be 
found  in  Appendix  B. 


2.2.3  Numerical  Integration.  In  this  research,  the  SGP4  propagator  will  be 
compared  to  a  numerical  integration  for  a  reference  satellite’s  orbit.  To  do  this,  the  above 
Hamiltonian  will  be  numerically  integrated  using  a  Hamming  fourth-order  predictor- 
corrector  algorithm  [7].  The  following  numerical  integrator  description  was  referenced 
from  Dr.  Wiesel’s  Hamming  Fortran  program.  First,  the  initial  conditions  of  the  orbit 
are  assigned  to  a  state  vector.  These  initial  conditions  are  the  generalized  coordinates 
and  momenta  at  the  start  of  the  integration  when  time  is  equal  to  zero, 


x  = 


=  < 


X 

y 

z 

X 

y 


> . 


(2.24) 
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Second,  the  derivative  of  the  state  vector  is  taken  and  the  Hamilton  canonical  equations 
are  substituted  in  for  q  and  p. 


x  = 


f  a  1 

f  P2L  ) 

J  q 

>  =  < 

dpk 

>  =  < 

l  p  J 

an 

^  9qk  J 

[ 

dV 

Py  ~  % 

dV 

U®Px  -  ^ 
_av 

dz 

Px  +  U®y 

Py  -  UJ^X 

Pz 


> . 


(2.25) 


Next,  the  total  time  of  the  integration  is  divided  by  the  number  of  steps  taken  during 
the  integration  to  get  an  incremental  time  step  ( dt) .  Finally,  the  future  position  is  the 
old  position  plus  the  change  in  position  multiplied  by  the  time  step, 


xi+1  =  Xi  +  x  idt.  (2.26) 

Numerical  integration,  though  accurate,  is  very  computationally  expensive  and  thus 
requires  large  amounts  of  computer  resources.  All  publicly  available  orbit  data  from  the 
DOD  Space-Track  website  is  propagated  with  SGP4.  However,  numerical  integration 
may  be  used  within  the  Air  Force  for  high  priority  objects,  collision  avoidance  situations, 
and  the  two  weeks  before  large  objects  re-enter  Earth’s  atmosphere  [4]. 

2.2.4  SGP4-  Simplified  General  Perturbation  4  (SGP4)  is  the  current  orbit 
propagator  employed  by  the  United  States  Air  Force.  It  was  developed  in  the  1960’s, 
went  operational  in  the  1970’s,  and  has  been  updated  many  times  during  its  history  [11], 
In  1988,  North  American  Aerospace  Defense  Command  (NORAD)  decided  to  release 
the  SGP4  code  giving  end  users  of  NORAD  two- line  element  (TLE)  sets  the  ability  to 
propagate  Earth  orbiting  objects  [12].  Since  that  time,  both  military  and  civilian  users 
updated  their  copy  of  SGP4  to  fit  a  specific  use.  This  created  many  working  versions  of 
the  code.  In  2006,  a  group  of  four  individuals  consolidated  these  working  versions  into 
a  standard  version  of  SGP4  that  hopefully  matched  the  Department  of  Defense  (DoD) 
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version  [11],  The  Air  Force  has  not  released  information  on  object  orbit  propagation 
since  1988,  so  it  is  impossible  to  know  for  sure  what  is  being  used. 

SGP4  treats  orbiting  objects  as  either  near-Earth  or  deep-space,  with  the  boundary 
at  an  orbital  period  of  225  minutes.  Earth’s  geopotential,  air  drag,  solar  pressure,  and 
lunar  effects  are  all  accounted  for  within  SGP4.  However,  the  1980  version  of  SGP4  used 
in  this  research  only  accounts  for  the  Earth’s  geopotential  to  order  and  degree  20  [11]. 

2. 3  Transformation  Theory 

2.3.1  Generating  Functions.  Hamiltonian’s  have  a  property  in  that  a  coordi¬ 
nate  transformation  can  be  made  to  a  different  Hamiltonian  with  different  coordinates 
and  momenta  while  still  describing  the  same  system.  These  transformations  provide 
an  opportunity  to  look  at  a  single  problem  multiple  ways.  Quite  frequently,  a  different 
Hamiltonian  will  have  coordinates  and  momenta  such  that  the  problem  is  solvable  or 
simplified  [13].  A  similar  idea  is  carried  out  when  a  person  in  a  room  is  identified  by 
saying  “look  at  two  o’clock”  instead  of  saying  “look  6  feet  east  and  8  foot  north” .  Both 
commands  point  in  relatively  the  same  direction,  but  one  gets  there  faster  and  easier. 

It  is  possible  to  switch  to  a  new  Hamiltonian,  but  how  and  which  one?  In  letters 
from  1836  and  1837,  Carl  Gustav  Jacobi  answered  this  question  with  the  Hamilton- 
Jacobi  transformation  theory  [14].  He  was  able  to  learn,  use,  and  improve  upon  the 
ideas  Hamilton  had  only  released  a  few  years  earlier. 

There  should  be  some  function,  F,  that  will  take  the  old  coordinates  and  momenta 
and  transform  them  into  the  new  coordinates  and  momenta.  At  this  point  neither  the 
function,  the  new  coordinates,  or  the  new  momenta  are  known,  but  the  statement  is  still 
accurate.  Therefore,  the  new  coordinates,  Q,  and  momenta,  P,  should  be  functions  of  the 
old  coordinates,  q,  and  momenta,  p,  and  the  new  Hamiltonian,  /C,  should  be  a  function 
of  the  new  coordinates  and  momenta  [13], 


Q 

Q(QmPnit)i 

(2.27) 

p 

=  P(qn,pn,t), 

(2.28) 

K 

=  IC(Q,P,t). 

(2.29) 
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With  one  set  of  old  and  one  set  of  new  coordinates  and  momenta,  there  are  four 
possible  combinations.  These  four  combinations  translate  to  four  different  functions  that 
transform  from  the  old  to  the  new  Hamiltonian.  Each  function  has  one  old  and  one  new 
coordinate  and  momenta,  and  with  this  information  the  remaining  two  coordinates  and 
momenta  can  be  determined  [13]. 

For  this  problem,  the  F2  generating  function  will  be  used  which  also  happens  to 
be  the  most  frequently  used  function.  The  properties  of  the  F2  generating  function  are 
listed  below  [15]. 


F2  = 

-^2 (?n j  Pn ?  5 

(2.30) 

dF2 

(2.31) 

Pi  = 

dqi  ’ 

d  F2 

(2.32) 

Qi  = 

dPl ' 

By  changing  the  coordinate  and  momenta  combinations  in  Equation  2.30,  the  other  three 
functions  and  relations  can  be  identified. 

At  this  point,  the  author  will  skip  forward  to  the  relationship  between  the  old 
and  new  Hamiltonians.  If  more  information  is  desired  on  generating  functions  or  the 
transformation  itself,  please  see  reference  [13],  pages  43-58.  Essentially,  the  derivative  of 
the  F-2  function  with  respect  to  time  is  broken  down  into  a  sum  of  its  partial  derivatives 
with  respect  to  the  old  coordinates,  separated,  and  substituted  in  for  the  new  momenta 
which  creates  the  generating  function.  This  also  explains  why  the  partial  of  the  F2  with 
respect  to  time  remains  [13]. 

2.3.2  Hamilton- Jacobi  Equation.  There  is  now  enough  information  to  form  a 
relationship  between  the  new  and  old  Hamiltonians  [13], 

OF 

fc(Qk,Pk)  =  ' H(qk,Pk )  +  ~wr-  (2.33) 
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Recall  that  when  a  coordinate  or  momenta  is  missing  from  the  Lagrangian  or  Hamilto¬ 
nian,  the  conjugate  momenta  or  coordinate  is  a  constant.  Therefore,  if  the  new  coordinate 
and  momenta  do  not  appear  in  the  new  Hamiltonian,  then  the  new  Hamilton  canonical 
equations  would  be  equal  to  zero  [13]. 


Qk 

die 

dPu 

(2.34) 

Pk  =  ~ 

dIC 

dQk 

(2.35) 

If  the  Hamilton  canonical  equations  are  zero  and  the  new  coordinates  and  momenta  are 
constant,  the  new  Hamiltonian  must  equal  zero.  This  changes  Equation  2.33  into  the 
Hamilton- Jacobi  equation  [13], 

OF 

'H(qk,Pk)  +  =  0.  (2.36) 

This  process  exentially  seeks  out  a  generating  function  such  that  the  new  coordinates 
and  momenta  are  integrals  of  the  motion  and  the  new  Hamiltonian  is  zero.  As  can  be 
expected,  this  transformation  greatly  simplifies  the  problem. 

2.4  KAM  Theory 

KAM  theory’s  beginnings  were  initially  proposed  by  A.N.  Kolmogorov  in  his  1954 
paper,  “On  conservation  of  conditionally  periodic  motions  for  a  small  change  in  Hamil¬ 
ton’s  function.”  In  the  1960’s,  J.  Moser  and  V.I.  Arnol’d  proved  the  theory  after  releasing 
multiple  papers  [9].  Relatively  speaking,  this  powerful  theory  has  gone  unused. 

For  a  completely  solvable,  or  integrable,  Hamiltonian,  all  solutions  lie  on  a  n- 
dimensional  torus  in  2 n- dimensional  phase  space,  where  n  is  the  number  of  coordinates. 
KAM  theory  simply  states,  lightly  perturbed  Hamiltonian  systems  will  have  most  of 
their  solutions  still  lying  on  the  surface  of  an  invariant  torus  given  the  perturbation  is 
small  enough  [16], [17], 

n  =H0(I)  +  en1{I,(j>),  e<l.  (2.37) 
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In  the  above  equation,  "H0  is  the  Hamiltonian  of  the  unperturbed  problem,  TL\  is  the 
perturbed,  periodic  Hamiltonian,  and  e  indicates  the  perturbations  are  very  small.  The 
above  Hamiltonians  for  the  KAM  torus  are  in  action-angle  variable  form,  meaning  action 
integral  variables,  I,  are  the  coordinates,  q’s,  and  angle  variables,  0,  are  the  conjugate 
momenta,  p’s.  Since  the  coordinates  and  momenta  are  different  than  our  current  Earth 
satellite  Hamiltonian,  Equation  2.19,  a  transformation  must  be  used.  The  process  dis¬ 
cussed  in  Section  2.3  concerning  generating  functions  and  the  Hamilton- Jacobi  equation 
will  be  used  and  discussed  more  in  the  following  section.  Notice  in  Equation  2.37  that  the 
unperturbed  system,  7/0(1),  is  missing  the  angle  variable  momenta,  <fi.  From  Equation 
2.12,  this  means  the  rate  at  which  the  coordinate  action  integral  variable,  I,  changes  is 
zero,  so  the  action  inegral  variable  is  a  constant. 

Invariant,  mentioned  above,  means  once  a  trajectory  in  phase  space  is  on  an  in¬ 
variant  torus  it  will  remain  on  the  torus.  As  an  example,  think  of  a  rolling  coin.  If  a 
coin  is  rolled  on  a  flat  surface  there  is  no  telling  which  direction  it  will  go  or  which  side 
it  will  fall  onto.  This  is  non- invariant.  Now,  think  of  the  spiraling  coin  funnels  found  in 
malls  and  museums.  Once  the  coin  is  release,  barring  some  outside  force  such  as  a  child 
picking  the  coin  up,  the  coin  will  continue  to  roll  around  in  circles  and  eventually  fall 
down  the  whole.  There  is  no  stopping  it  and  this  is  invariant.  The  only  difference  is  the 
coin  eventually  stops,  but  the  satellite  motion  on  the  KAM  torus  is  periodic  so  it  will 
continue  repeating  forever.  If  an  Earth  satellite  is  to  lie  on  a  KAM  torus,  it  will  do  so 
in  phase  space  and  it  will  be  a  lightly  perturbed  invariant  torus  [7]. 

2-4-1  Transformation.  Discuss  here  after  talking  with  Wiesel 

Earlier,  transformation  theory  and  the  Hamiltonian- Jacobi  equation  were  discussed. 
Equation  2.19  is  the  Hamiltonian  for  an  Earth  orbiting  satellite  with  coordinates  x,  y, 
and  z.  The  KAM  theory  Hamiltonian,  Equation  2.37,  is  in  action  and  angle  variables. 
To  see  if  Earth  orbits  can  be  described  using  KAM  tori,  it  is  necessary  to  transform  the 
Earth  satellite  Hamiltonian,  Equation  2.19,  into  the  form  of  the  KAM  theory  Hamilto¬ 
nian,  Equation  2.37.  This  can  be  done  using  the  Hamilton- Jacobi  equation,  2.36,  and 
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methods  discussed  in  Section  2.3  on  transformation  theory.  Dr.  Wiesel  has  shown  this 
process  in  reference  [16]. 

2-4-2  KAM  Torus  Orbital  Coordinates.  Dr.  Wiesel  applied  a  Fast  Fourier 
Transform  (FFT)  algorithm  to  a  one  year  numerical  integration  of  the  reference  satellite. 
The  algorithm  was  able  to  separate  the  periodic  behavior  of  the  orbit  into  three  base 
frequencies.  Dr.  Wiesel  shows  in  [16]  that  these  three  frequencies  used  to  describe  the 
KAM  torus  can  be  linked  to  three  classical  Delaunay  orbital  elements.  The  largest  of 
these  three  frequencies  is  the  ‘mean’  mean  anomaly,  M.  The  second  frequency  is  the 
longitude  of  the  mean  node,  Q.  The  smallest  frequency  is  the  mean  argument  of  perigee, 
uj.  These  three  angles,  with  the  use  of  Hamilton- Jacobi  theory,  will  be  transformed  from 
the  x,  y,  z  positional  coordinates  of  the  original  Hamiltonian  [7]  to  the  coordinates  of  the 
new  KAM  Hamiltonian. 

2-4-3  Phase  Space  KAM  Torus.  In  phase  space,  the  three  fundamental  fre¬ 
quencies  listed  above  creates  a  torus  of  a  particular  size  and  shape  within  2n-dimensional 
phase  space  that  describes  the  satellite’s  orbit.  The  torus  is  described  with  a  Fourier 
series  including  these  frequencies.  Also,  the  coordinates  and  momenta  are  grouped  to¬ 
gether  into  one  state  vector.  During  the  numerical  integration  of  the  satellites  orbit,  the 
rate  of  change  of  the  state  vector  is  found.  Using  a  least  squares  fitting  algorithm,  the 
state  vector  and  its  rate  of  change  can  be  made  to  fit  a  Fourier  series  torus  that  includes 
the  three  fundamental  frequencies  describing  the  orbit.  Once  converged,  this  gives  all 
the  coefficients  necessary  to  uniquely  describe  the  satellite  as  a  Fourier  series,  which  in 
turn  produces  a  torus  in  2n-dimensional  phase  space  [7]. 
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III.  Methodology 

This  section  lays  out  the  processes  taken  during  the  course  of  this  research.  Since  the 
research  was  computer  simulation,  computer  and  software  setup  will  be  discussed  first. 
Any  assumptions  made  during  this  research  will  then  be  listed.  Finally,  the  programing 
and  testing  process  will  be  discussed  for  SGP4,  STK,  and  KAM  torus  orbit  propagation 
techniques. 

3.1  Computer  Software 

3.1.1  Fortran  Programming.  The  author  performed  all  programming  within 
the  Microsoft®  Visual  Studio®  2008  program  development  environment  on  a  Windows 
Vista®  64-bit  computer.  The  compiler  was  Intel®  Visual  Fortran  Compiler  Professional 
Edition  for  Windows®,  version  11.1.046.  Both  software  suites  were  used  with  educational 
licenses.  Almost  all  Fortran  code  was  previously  written  and  provided  to  the  author  by 
Dr.  Wiesel.  Code  was  updated  to  Fortran  95/2003  wherever  possible  and  edited  to 
run  within  the  Microsoft®  Visual  Studio®  (MVS)  environment.  Double  precision  was 
chosen  for  all  real  variables  with  a  user  defined  length  of  15  digits.  This  helps  ensure 
maximum  portability  and  consistency  on  different  computer  systems  and  allows  for  future 
precision  changes  from  one  location.  The  ‘IMPLICT  NONE’  command  was  always  used 
to  minimize  possible  variable  errors. 

3.1.2  Commercial  Orbit  Prediction.  In  order  to  help  validate  previous  results, 
it  was  necessary  to  compare  to  an  outside  source.  Analytical  Graphics®  (AGI)  Satellite 
Toolkit®  (STK)  was  easily  accessible  with  a  free  educational  license  and  was  recommend 
by  Dr.  Wiesel.  The  author  used  STK  Version  9.0.1,  which  has  multiple  propagators 
that  will  determine  a  satellite’s  position  in  time.  The  author  choose  the  High  Precision 
Orbit  Propagator  (HPOP)  due  to  its  ease  of  use  and  detailed  capabilities  far  exceeding 
necessary  requirements.  The  HPOP  propagator  is  a  numerical  integrator  utilizing  the 
differential  equations  of  motion  for  the  satellite  initially  generated  from  ephemeris  data, 
as  stated  in  the  STK  help  hies.  For  those  unfamiliar,  ephemeris  data  includes  position 
and  velocity  of  a  satellite  for  a  given  time  period.  STK  allows  the  user  to  import  a  specific 
satellite  hie,  choose  the  date  and  time,  and  ephemeris  data  is  automatically  generated 
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via  an  AGI  server  connection.  Almost  any  force  or  perturbation  can  be  modeled  with 
the  HPOP  propagator  in  great  detail. 

3.1.3  Plotting.  All  plotting  was  completed  using  the  Windows®  version  of 
Gnuplot©,  Version  4.2.6.  The  software  is  free  and  relatively  easy  to  use  for  simple 
graphs.  However,  the  program  has  many  advanced  features  allowing  for  very  detailed 
and  specific  plots.  Dr.  Wiesel  uses  Gnuplot  and  example  scripts  were  easily  edited  by 
the  author. 

3.2  Assumptions 

Given  the  academic  nature  and  initial  testing  of  this  theory,  certain  effects  and 
details  have  not  been  accounted  for  that  normally  would  for  a  developed  technique  used 
by  NASA  on  an  actual  satellite.  The  geopotential,  Equation  2.15,  discussed  in  Section 
2.2.1  is  used  to  calculate  the  Earth’s  gravitational  force  on  the  satellite  as  it  orbits  the 
Earth.  For  this  research,  order  and  degree  truncated  to  20  was  decided  sufficient. 

Air  drag,  ocean  tides,  lunar  and  solar  gravitational  effects,  and  other  third  body 
effects  have  all  been  assumed  small  enough  when  compared  to  Earth’s  gravity  to  be  ne¬ 
glected.  If  these  smaller  terms  had  been  accounted  for  then  a  more  detailed  geopotential 
model,  higher  order  and  degree  for  Equation  2.15,  would  have  been  used. 

While  actual  satellite  mission  orbit  prediction  would  use  most  if  not  all  the  above 
effects,  the  mission  of  this  research  to  show  KAM  theory  can  be  applied  to  Earth  satellites 
for  orbit  prediction  can  reliably  succeed  without  these  effects. 

3.3  Reference  Satellite 

Continuing  with  Dr.  Wiesel’s  work,  the  reference  satellite  for  this  research  will  be 
the  same.  This  will  allow  for  continuity  between  papers  and  with  future  research.  The 
reference  satellite’s  initial  conditions  are  listed  in  Table  3.1,  orbital  elements  are  listed 
in  Table  3.2,  and  a  picture  of  the  orbit  is  in  Figure  3.1. 

3.4  SGP4  Verification 
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Table  3.1:  Reference  Satellite  Initial  Conditions 


X 

y 

z 

Position  (km) 

-4412.83115168178 

4676.00408732872 

-2910.15168627727 

Velocity  (km/sec) 

-4.42110219123704 

-5.02218008450107 

-2.27614671171868 

Table  3.2:  Reference  Satellite  Orbital  Elements 


Semi- major 
Axis 

Eccentricity 

Inclination 

Right  Ascension 
Ascending  Node 

Argument  of 
Perigee 

Mean 

Anomaly 

a 

e 

i 

O 

id 

M 

7,049.5  km 

0.05 

30.0° 

261.72° 

141.41° 

88.42° 

Figure  3.1:  Reference  Satellite  Orbit  from  STK. 


3-4-1  Test  Case.  The  version  of  SGP4  used  is  in  this  research  is  different  from 
the  original  1980  version  of  SGP4  released  by  Hoots  and  Roehrich  in  [12].  As  stated 
earlier,  the  author  was  given  Dr.  Wiesel’s  existing  Fortran  code.  The  first  step  was  to 
make  sure  it  would  run  properly  on  the  author’s  computer  and  produce  correct  results. 
Starting  with  the  SGP4  program,  a  test  satellite  was  both  numerically  integrated  using 
a  Hamming  fourth-order  predictor-corrector  algorithm  and  propagated  with  SGP4.  The 
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differences  in  x,  y,  z  position  between  the  numerical  integrator  and  SGP4  results  are 
called  the  residuals.  These  residuals  are  then  plotted  both  individually  and  together  as 
a  total  magnitude. 

3.Jh2  Reference  Satellite.  For  comparison  later  in  this  paper,  the  above  process 
was  also  completed  for  the  reference  satellite  discussed  above.  Throughout  this  research 
the  reference  satellite  was  used  to  help  ensure  proper  program  operation  and  for  result 
comparison  with  other  satellites  that  were  used. 

3.5  Comparing  SGP4  and  STK  Propagators 

Next,  the  SGP4  satellite  propagator  was  compared  to  STK’s  HPOP  propagator. 
The  satellites  tested  were  randomly  chosen  from  all  orbiting  satellites  with  only  two 
requirements.  The  first  being  a  relatively  equal  spacing  of  satellites  in  the  low  Earth 
orbit  (LEO)  region.  Second,  within  the  LEO  group  a  few  satellites  had  to  be  very  close 
to  each  other  to  compare  results.  Satellites  from  medium  Earth  orbit  (MEO)  and  high 
Earth  orbit  (HEO)  were  not  included  in  this  test  since  the  deep-space  section  of  the 
SGP4  code  was  not  included  in  the  version  used  in  this  research.  The  satellites  chosen 
are  listed  in  the  following  table  with  data  taken  from  STK  at  the  initial  time  prior  to 
propagating.  The  satellites  used  in  the  SGP4  versus  STK  comparison  can  be  seen  in 


Table  3.3:  STK  Satellite  Orbital  Data  &  Elements 


Satellite 

Altitude 

(km) 

Semi-Major  Axis 
(km) 

Eccentricity 

Inclination 

(degrees) 

Period 

(mins) 

Grace  1 

467.94 

6839.52 

0.00125 

88.99° 

93.82 

Grace  2 

468.23 

6838.84 

0.00133 

88.99° 

93.81 

Swift 

584.6 

6972.19 

0.00196 

20.53° 

96.56 

Reference  Satellite 

679.28 

7049.5 

0.05 

29.96° 

98.17 

Figure  3.2. 

Each  satellite  had  position  and  velocity  information  created  with  HPOP  which 
was  read  into  the  SGP4  and  STK  Fortran  program.  The  SGP4  propagator  used  the  first 
line  of  each  satellite’s  inertial  data  as  the  initial  conditions  and  proceeded  from  there. 
All  tested  satellites  were  initiated  with  the  same  start  date  and  time.  Satellites  were 
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Figure  3.2:  Satellites  used  in  SGP4  /  STK  comparison. 


propagated  out  in  increments  of  10,  15,  30  and  60  days.  The  residuals,  the  difference  in 
location  between  the  two  propagators,  were  calculated  at  each  time  step.  Each  satellite 
that  completed  the  program  and  showed  residuals  of  tens  of  meters  or  less  would  then 
be  run  with  the  next  higher  propagation  time  step. 

3.6  Torus  Frequencies 

An  Earth  satellite  orbit  can  be  approximated  as  a  KAM  torus  [7].  The  primary 
purpose  of  this  research  is  to  verify  this  statement  against  simulated  STK  orbit  data. 
Dr.  Wiesel  numerically  integrated  a  reference  satellite  for  one  year,  with  x,  y,  z  position 
data  stored  in  a  hie  for  each  time  step.  Reference  satellite  characteristics  were  discussed 
in  Section  3.3,  and  the  process  for  Ending  the  three  torus  frequencies  was  discussed  in 
Section  2.4.2. 
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A  program  written  by  Dr.  Wiesel  called  ‘PeakEater’  takes  an  orbit  file  with  position 
and  time  data  and  produces  a  file  describing  a  torus.  The  torus  is  created  using  a  Fourier 
series, 

n 

cos (ji  ■  Q)  T  Sj  sin(ji  ■  Q),  (3.1) 

3= 1 

where  Cj  and  Sj  are  the  Fourier  series  coefficient  vectors,  the  Q  vectors  is  the  three  torus 
frequencies  multiplied  by  time  (Qi  =  c 0{t  +  Qoi),  and  j  is  an  indice  vector  indicating 
each  torus  frequency  has  its  own  set  of  Fourier  series  coefficients.  Remember,  the  torus 
has  n-dimensions  in  2n-dimensional  phase  space. 

Once  the  location  on  a  torus  is  found,  that  point  linearly  drifts  with  time, 


Q  =  { 


Uit  +  Q  oi 
U}2t  +  Q  02 
^3 1  +  Qo3 


(3.2) 


This  linear  drift  around  the  torus  in  phase  space  is  the  satellite  propagating  around  the 
orbit  in  real  space. 

Using  the  reference  satellite’s  initial  conditions,  STK  produced  a  120  day  orbit 
propagation  with  HPOP.  A  program  from  Dr.  Wiesel  called  ‘Coordinate  Residuals’ 
compared  the  STK  orbit  to  the  KAM  torus  Fourier  series.  Once  the  location  of  the  STK 
orbit  is  found  on  the  torus,  the  orbit  is  propagated  forward  using  the  time.  The  location 
on  the  torus  is  then  converted  back  into  positional  coordinates  and  these  are  compared 
to  the  STK  position  to  again  produce  residuals.  These  residuals  are  again  plotted  as 
individual  coordinate  residuals  and  as  a  residual  magnitude. 
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IV.  Results  and  Discussion 


In  this  section  the  results  will  be  presented  and  discussed.  As  a  reminder,  coordinate 


residuals  are  the  differences  between  each  coordinate  (x,  y,  z)  for  the  two  methods  being 


compared.  Residual  magnitude  is  the  magnitude  of  the  residual  coordinates  for  the  two 


methods  being  compared. 

4-1  SGP4  Verification 

4-1.1  Test  Case.  The  first  step  was  to  ensure  all  code  from  Dr.  Wiesel  was 
implemented  correctly.  After  running  the  test  case,  output  plots  of  coordinate  residuals 
and  residual  magnitude  of  all  three  coordinates  were  compared.  The  initial  plots  using 
SGP4  matched  Dr.  Wiesel’s  and  thus  confirmed  the  initial  SGP4  propagator  and  nu¬ 
merical  integrators  were  working.  Figure  4.1  shows  the  x,  y,  and  z  residuals  between  the 
numerical  integration  and  the  SGP4  propagator  stayed  under  1.2  kilometer  for  the  first 
6  days,  while  Figure  4.2  shows  the  magnitude  of  the  residual  error  remained  under  1.25 
kilometers  over  the  6  day  period. 
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Figure  4.1:  SGP4  verification  showing  coordinate  residuals  for  test  case. 
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Figure  4.2:  SGP4  verification  showing  residual  magnitude  for  test  case. 

4-1.2  Reference  Satellite.  With  the  SGP4  program  working  properly,  the  ref¬ 
erence  satellite  was  run  next  for  a  fifty-five  day  period.  Figure  4.3  shows  individual 
coordinate  residuals  as  high  as  six  kilometers  and  Figure  4.4  shows  error  magnitude  at 
a  maximum  of  seven  kilometers  during  the  55  day  period.  These  plots  will  be  used  for 
comparison  later  in  the  paper. 

4-2  Comparing  SGP4  and,  STK  Propagators 

The  results  of  comparing  the  SGP4  and  STK  propagators  are  discussed  here.  The 
Grace  1  and  Grace  2  satellites  were  chosen  for  their  similar  orbits.  Results  for  these  two 


Table  4.1:  SGP4/STK  Satellite  Residual  Results 


Satellite 

Comparison 

Length 

(days) 

Coordinate  Residuals 

Residual 

Magnitude 

(km) 

X 

(km) 

y 

(km) 

z 

(km) 

Grace  1 

25 

3.5 

4 

3 

4.5 

Grace2 

20 

3 

3.5 

2.5 

4 

SWIFT 

25 

2 

2 

3 

3 

Reference  Satellite 

55 

3 

2 

3 

3.5 
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Figure  4.3:  SGP4  verification  showing  coordinate  residuals  for  reference  satellite. 


Figure  4.4:  SGP4  verification  showing  residual  magnitude  for  reference  satellite. 

satellites  should  be  very  similar.  As  can  be  seen  in  Figures  4.5  to  4.8,  the  results  are 
indeed  very  similar  and  help  to  show  the  program  is  working  correctly.  If  results  for 
these  two  satellites  had  been  off,  the  author  would  have  been  alerted  to  look  for  an  error 
that  otherwise  may  have  gone  unnoticed  if  all  satellites  were  dissimilar. 
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Figure  4.5:  SGP4/STK  coordinate  residuals  for  Grace  1  satellite. 
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Figure  4.6:  SGP4/STK  residual  magnitude  for  Grace  1  satellite. 

The  comparison  between  SGP4  and  STK  for  the  Grace  1  and  Grace  2  satellites 
was  good  to  under  4.5  and  4  kilometers  in  magnitude  respectively.  A  complete  listing  of 
coordinate  residuals  and  residual  magnitudes  can  be  seen  in  Table  4.1  and  plots  for  the 
remaining  satellites  can  be  found  in  Appendix  D. 
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Figure  4.7:  SGP4/STK  coordinate  residuals  for  Grace  2  satellite. 


Figure  4.8:  SGP4/STK  residual  magnitude  for  Grace  2  satellite. 

Some  satellites  tested  would  not  produce  residuals  under  100  kilometers  within 
the  first  ten  to  fifteen  days.  The  actual  reason  is  not  know  at  this  time,  however  there 
are  some  possibilities.  The  original  1980  SGP4  code  did  have  some  errors  for  certain 
satellite  orbits.  This  was  one  of  the  reasons  users  started  modifying  the  1980  version. 
Another  possibility  is  the  chance  of  errors  within  the  least  squares  converging  process. 
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The  matrices  can  be  singular  or  near  singular,  with  high  condition  numbers,  which  could 
cause  problems  in  convergence. 

4-3  STK  and  Torus  Comparison  for  Reference  Satellite 

As  mentioned  earlier,  the  reference  satellite  was  numerically  integrated  by  Dr. 
Wiesel  for  one  year.  This  one  year  long  numerical  integration  was  used  to  identify  the 
three  torus  frequencies  which  describe  the  reference  satellite  as  a  torus  in  phase  space. 
The  reference  satellite  was  propagated  for  60  days  with  STK.  This  60  days  of  orbit  data 
was  then  fit  to  the  year  long  reference  torus.  The  differences  in  individual  coordinates 
and  total  magnitude  were  calculated  as  before  and  used  for  plotting. 

Initial  residuals  reached  500  kilometers  within  60  days  at  a  linear  rate  as  seen  in 
Figure  4.9. 
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Figure  4.9:  Torus  vs.  STK  data  coordinate  residuals  for  reference  satellite. 


Dr.  Wiesel  suggested  the  frequencies  describing  the  numerically  integrated  orbit 
may  differ  from  the  STK  orbit  frequencies  since  their  calculation  methods  are  different. 
If  the  frequencies  used  to  fit  the  STK  orbit  data  to  the  reference  satellite  torus  were 
adjusted,  better  accuracy  could  be  achieved.  Of  the  three  torus  frequencies,  ay,  the 
mean  motion,  is  the  largest  and  would  be  the  one  to  adjust.  First,  the  slope  with  which 
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Figure  4.10:  Torus  vs.  STK  data  residual  magnitude  for  reference  satellite. 

the  residual  magnitude  rose  in  Figure  4.10  was  found  at  multiple  points  and  had  units 
of  kilometers  per  day.  This  was  done  to  get  an  average  value  for  the  entire  plot.  Each 
calculated  slope  was  then  divided  by  the  distance  the  point  used  was  from  the  origin 
of  the  ECEF  system.  This  then  gave  the  amount  of  radians  per  day  the  residuals  were 
increasing.  Finally,  days  were  converted  to  time  units  which  are  used  within  the  program. 
The  change  in  frequency  was  rather  small  at  a  0.0012%  change.  The  average  value  for 
the  plot  was  the  frequency  difference  between  the  numerically  integrated  torus  and  the 
STK  orbit  data  torus.  By  trial  and  error  it  was  determined  if  this  value  should  be  added 
or  subtracted  from  the  original  torus  frequency, 

After  implementing  the  frequency  change  and  re-running  the  programs,  residuals 
improved  to  thirty  kilometers  within  sixty  days,  or  an  over  90%  improvement.  The  entire 
process  was  completed  a  second  time  with  the  frequency  only  changing  by  6  X  10-5%. 
However,  the  residuals  again  improved  to  under  5  kilometers  in  a  60  day  period  which 
was  another  improvement  of  80%.  These  results  can  be  seen  in  Figures  4.11  to  4.12. 

Finally,  STK  orbit  data  for  the  reference  satellite  was  propagated  120  days  and  fit 
to  the  reference  satellite  torus.  As  seen  in  the  following  plots,  coordinate  residuals  were 
under  three  kilometers  and  residual  magnitude  remained  under  3.5  kilometers. 
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Figure  4.11:  Torus  vs.  STK  data  coordinate  residuals  for  reference  satellite. 
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Figure  4.12:  Torus  vs.  STK  data  residual  magnitude  for  reference  satellite. 

Recalling  Section  4.2,  SGP4  was  able  to  compare  to  STK  orbit  data  for  the  refer¬ 
ence  satellite  out  to  55  days  with  3.5  kilometers  of  residual  magnitude  error,  but  could 
not  converge  on  anything  greater  than  sixty  days.  The  torus  model  allows  much  great 
comparison  time  periods  with  equal  or  better  accuracy. 
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Figure  4.13:  Torus  vs.  STK  data  coordinate  residuals  for  reference  satellite. 


Figure  4.14:  Torus  vs.  STK  data  residual  magnitude  for  reference  satellite. 

It  is  also  important  to  note  that  when  working  in  the  ‘fixed’  frame  in  STK,  this 
paper’s  ECEF  frame,  some  perturbations  not  included  in  this  paper  have  been  added 
back  in  automatically.  STK  is  not  intended  to  be  used  as  a  simplified  model,  so  some 
perturbations  were  ‘hard-coded’  into  STK.  These  ‘hard-coded’  perturbations  could  be 
accounting  for  the  patterns  seen  in  the  figures  presented  in  this  chapter. 
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V.  Conclusions 


5.1  Comparing  SGP4  and  STK 

The  results  from  comparing  actual  satellites  with  SGP4  and  STK  propagation 
methods  produce  relatively  low  residuals.  In  normal  location  monitoring,  times  between 
tracks  would  probably  not  reach  twenty  to  twenty-five  days  as  seen  here.  Therefore  the 
three  to  five  kilometers  in  magnitude  error  would  be  plenty  low  enough  to  ensure  that 
satellite  was  not  lost.  These  results  also  show  how  quickly  predicted  results  can  become 
large  and  how  frequently  these  predictions  must  be  updated  to  ensure  satellite  location 
is  known  within  a  reasonable  distance  and  reliability.  Applying  these  techniques  to  a 
17,000  object,  and  growing,  space  catalog  is  indeed  a  very  daunting  task. 

5.2  KAM  Torus  Orbit  Prediction 

The  results  obtained  from  comparing  STK  orbit  data  with  a  known  KAM  torns 
model  were  surprisingly  accurate.  At  time  periods  out  to  120  days,  individual  coordinate 
residuals  remained  under  three  kilometers.  Time  periods  greater  than  120  days  were 
tried  and  could  have  been  presented  here.  However,  at  those  time  lengths  STK  orbit 
reliability  has  diminished  with  the  simplified  model  being  used  in  this  research  and  any 
results  would  not  have  as  much  meaning  as  the  more  accurate  120  day  orbit  periods. 

With  SGP4  and  STK  orbit  prediction  techniques,  errors  slowly  start  to  creep  into 
results  and  grow  as  orbit  propagation  time  periods  increase.  However,  the  KAM  torus 
model  showed  constant  coordinate  residuals.  This  shows  the  periodic  nature  of  the  KAM 
torus  and  the  linear  drift  of  the  position  of  the  satellite  on  the  torus  in  phase  space  keeping 
the  residuals  at  constant  low  values.  It  also  shows  the  theory  behind  this  KAM  torus, 
that  lightly  perturbed  Hamiltonian  system’s  solutions  he  on  a  KAM  torus  and  once  on 
this  invariant  torus  the  solution  stays  on  it. 

5. 3  Recommendations 

As  part  of  the  initial  studies,  this  research  has  done  very  well.  The  KAM  torus 
model  takes  into  account  the  entire  geopotential  model  while  other  existing  techniques 
must  estimate  the  geopotential.  However,  the  KAM  torus  method  does  not  yet  include 
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other  perturbing  forces  that  are  commonly  accounted  for  in  orbit  prediction.  A  future 
step  would  be  to  incorporate  effects  such  as  air  drag,  lunar  effects  ,  and  solar  effects  and 
verify  if  the  solution  to  the  Hamiltonian  still  lies  on  the  surface  of  a  torus. 

It  is  believed  Dr.  Wiesel  will  be  looking  into  satellites  drifting  from  one  invariant 
torus  to  a  nearby  invariant  torus,  so  this  topic  holds  avenues  for  possible  research. 

Finally,  it  would  be  outstanding  to  work  with  an  existing  satellite  team  whose 
satellite  has  been  accurately  modeled  with  current  KAM  torus  methods  and  work  towards 
updating  the  KAM  torus  method  to  the  degree  that  the  satellite  team  would  navigate 
their  satellite  with  the  KAM  torus  method. 

Also,  Dr.  Wiesel  has  discussed  using  KAM  torus  methods  to  navigate  a  constella¬ 
tion  of  satellites.  With  the  linear  drift  of  a  satellite  on  the  surface  of  a  torus  in  phase 
space,  station  keeping  between  satellites  would  merely  entail  maintaining  equal  linear 
spacing  on  the  surface  of  the  torus  in  phase  space. 

Although  KAM  theory  is  relatively  old,  KAM  torus  orbit  prediction  methods  are 
relatively  new  and  hold  enormous  possibilities  for  the  future.  Hopefully  these  ideas  are 
getting  more  people  interested  and  the  rate  of  research  and  testing  can  increase. 
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Appendix  A.  Hamiltonian  Mechanics  Example 

In  this  appendix,  the  simple  pendulum  will  be  used  as  an  example  to  illustrate  the 
application  of  Lagrangian  and  Hamiltonian  mechanics  discussed  in  Chapter  II.  A  simple 
pendulum  is  fixed  to  an  upper,  non-moving  surface.  The  pendulum  is  allowed  to  move 
sideways  in  the  plane  of  the  paper  about  the  point  where  it  meets  the  surface. 


Figure  A.l:  A  Simple  Pendulum 


As  discussed  earlier,  multiple  coordinate  systems  can  usually  be  used  for  a  given 
problem.  For  this  problem  Cartesian  coordinates  can  be  used,  but  polar  coordinates 
simplify  the  problem, 


Q  = 


l 

9 


(A.l) 


The  kinetic  energy  of  the  pendulum  is, 


(A.2) 


where  the  velocity  for  the  pendulum  is  the  angular  velocity,  v  =  16.  This  results  in  the 
kinetic  energy  being, 

T  =  ]anl292.  (A.3) 

The  potential  energy  of  the  pendulum  based  on  its  height  is, 


V  =  mgh. 


(A.4) 
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The  height  of  the  pendulum  mass  can  be  determined  from  its  vertical  distance,  h  = 
—l  cos  9,  which  gives  a  potential  energy  of, 

V  =  mg(— l  cos  9)  =  —mglcos9.  (A. 5) 


As  mentioned  in  Section  2.1.3,  the  Lagrangian  of  a  system  is  the  kinetic  energy 
minus  the  potential  energy,  and  for  this  system  is, 

£  =  T  —  V  =  ml20 2  —  (— mgl  cos  9)  =  ^ ml29 2  +  rngl  cos  6.  (A. 6) 

To  simplify  the  problem,  the  Lagrangian  per  unit  mass  will  be  used  instead, 

£  =  —  =  -l292  +  gl  cos  9.  (A. 7) 

m  2 

Remembering  the  generalized  velocities  are, 


Q  = 


i 

9 


(A.8) 


and  the  Lagrangian  equations  of  motion  (EOM)  come  from, 

d(d£\_d£ 
dt\dqk)  dqk 

the  Lagrangian  equations  of  motion  (EOM)  can  be  solved  for, 


(A.9) 


l2 9  +  gl  sin  6  =  0. 


(A.10) 


However,  in  this  example  the  Hamiltonian  is  being  solved.  The  generalized  momenta  is, 


The  Hamiltonian  comes  from, 

n 

U  =  ^Pk4k  -  £■ 
k=  1 

Solving  for  the  generalized  velocity  from  the  generalized  momenta, 


9  = 


Pe 

Z2’ 


and  substituting  into  the  Hamiltonian  gives, 


U  =  -  gl  cos  9). 


The  Hamilton  canonical  equations  come  from, 
dU 


*  =  (A'15) 


and  are, 


Pk  = 


dU 

dqk 


(A. 16) 


A_Pe 

e~¥’ 


(A.17) 


Pe  =  —gls'm9.  (A. 18) 


(A. 12) 


(A. 13) 


(A. 14) 
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Appendix  B.  Earth  Orbiting  Satellite  Hamiltonian  Derivation 

In  this  appendix,  the  Hamiltonian  for  the  Earth  orbiting  reference  satellite  used  through¬ 
out  this  paper  will  be  derived.  The  same  steps  seen  in  the  simple  pendulum  example  in 
Appendix  A  will  be  followed. 


For  the  Earth  orbiting  satellite  problem,  an  Earth- Centered,  Earth-Fixed  (ECEF) 
frame  of  reference  will  be  used.  Regular  Cartesian  coordinates  (x,  y,  z)  will  be  used  to 
describe  position, 

M 


<?  =  < 


y 


(b.i) 


l  *  J 


The  velocity  of  the  satellite  is  simply  the  time  derivative  of  the  position.  However, 
remember  the  ECEF  reference  frame  is  rotating,  so  to  get  the  true  velocity  an  additional 
term  must  be  added  to  take  into  account  this  rotation  about  the  z-axis. 


Q  = 


x  — 

<  y  +  uj^x 


(B.2) 


The  kinetic  energy  can  then  be  determined  by  simply  plugging  the  velocity  into  Equation 
2.4, 

T  =  ^m((x  -  uj^y)2  +  (y  +  taex)2  +  i2).  (B.3) 

For  simplicity,  it  is  not  uncommon  to  divide  the  above  equation  through  by  the  mass 
which  then  gives  kinetic  energy  per  unit  mass.  One  reason  for  doing  this  is  not  being 
required  to  carry  the  mass  throughout  then  entire  process  that  follows,  while  another  is 
the  mass  may  not  yet  me  know  for  the  satellite.  Either  way,  in  this  paper  the  kinetic 
energy  and  all  quantities  with  mass  that  follow  will  be  per  unit  mass. 


T 


w<$yY  +  (y  +  ^®XY  +  z2). 


(B.4) 
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The  potential  energy,  as  discussed  in  Section  2.2.1,  is  given  by, 


V=-^YJ  V'!,  )  +  Snm  sinmA).  (B.5) 

The  Lagrangian  is  then  found  by  plugging  the  kinetic  and  potential  energy  into  equation 
B.6, 


1  oo  n 

£  =  ^((x-u9y)2+(y+uj9x)2+z2)--  ^  5)(Cnm  cosm\+Snmsmm\). 

2  r  n=l  m=l 

(B.6) 


The  canonical  momenta  are  obtained  using  equation  2.8, 


(  > 

/ 

Px 

x-uj®y 

p  =  < 

Py 

>  =  < 

y  +  ^®x 

p* 

V.  > 

Z 

\ 

Finally,  the  Hamiltonian  is  obtained  using  equation  2.10, 


(B.7) 


n 


1 

2 


(Px  +  Py  +  P2:)2  +  U)f$(ypx 


XPy ) 


-  r  (^)  "-^(sin  5) (Cnm  cos mA  +  sin mA) 

n=l  ra=l 


(B.8) 
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Appendix  C.  Non-linear  Least  Squares  Estimation  within  SGP4 

C.l  Least  Squares  Origin 

The  least  squares  estimation  process  is  a  method  of  using  observations,  or  data 
points,  to  accurately  predict  the  next  data  point  for  a  system.  Consider  a  non-linear 
system  for  which  the  solution  is  a  function  of  time, 

x(t)  =  h(x(t0),t).  (C.l) 

If  the  initial  location  and  time  are  know,  the  current  location  can  be  found.  In  this 
paper,  the  above  system  is  the  orbiting  satellite  being  propagated  via  SGP4. 

The  data  the  author  would  like  SGP4  to  model  is  the  STK  orbit  data  which  will  be 
called  the  observations,  z.  For  each  step  through  time,  SGP4  results  will  be  compared  to 
the  STK  data.  The  difference  between  the  two  is  an  error,  since  the  ideal  case  would  be 
for  SGP4  to  exactly  match  STK  and  produce  no  error.  The  error  is  the  SGP4  location 
minus  the  STK  location  for  each  time  step, 

e  =  x  —  z.  (C.2) 

These  errors  are  then  added  up  to  determine  the  total  error  produced  during  the  prop¬ 
agation, 

N 

£  =  X>-  (C.3) 

i= 1 

The  case  in  which  SGP4  propagation  best  matched  the  STK  data  should  create  the 
smallest  total  error,  E.  However,  if  the  SGP4  data  was  smaller  than  the  STK  data  for  a 
given  time  step  a  negative  number  would  be  produced.  When  added  up  at  the  end,  this 
would  make  the  total  error  look  smaller  even  though  it  was  not.  To  fix  this  problem  the 
individual  errors  for  each  time  step  are  squared,  ensuring  all  numbers  are  positive, 

e2  =  (x  —  z)2.  (C.4) 
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These  squared  errors  are  then  added  up  at  the  end  to  produce  a  total  squared  error, 


N 

£2  =  £^'  (C.5) 

i= 1 

Even  though  the  errors  have  been  squared,  the  case  in  which  SGP4  and  STK  matched 
the  best  would  still  be  when  the  squared  error  was  the  least.  This  is  where  least  squares 
estimation  gets  its  name. 


C.2  Non-Linear  Least  Squares 


Starting  with  Equation  C.l,  SGP4  will  propagate  a  satellite  through  time  non- 
linearly.  As  mentioned  earlier,  this  is  completed  by  incrementing  the  orbital  elements 
through  time  to  produce  a  new  location  and  velocity.  The  state  of  the  satellite  will  be 
used  to  refer  to  its  position  and  velocity, 


x  =  < 


X 

y 

z 

X 

y 


(C.6) 


In  an  effort  to  predict  the  future  state  of  the  system,  the  dynamics  (Equation  C.l)  will 
be  linearized.  Essentially,  to  get  from  the  starting  point  to  the  ending  point  a  bunch 
of  small  steps  will  be  taken.  To  get  from  one  state  to  the  next  state,  a  state  transition 
matrix,  <E>,  is  use, 

Sx(ti+1)  =  ®(ti+1,ti)5x(ti),  (C.7) 


where  the  next  state  is, 


x(ti+ 1)  =  Xi  +  Sx(ti+ 1), 


(C.8) 


Earlier  it  was  said  the  SGP4  process  is  non-linear.  The  entire  SGP4  process  will  be 
represented  by  a  non-linear  function,  G,  such  that  the  state  (Equation  C.6)  and  time 
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put  through  this  function  will  produce  the  correct  STK  orbit  data  (z)  for  that  time, 


Zi(ti)  =  G{x{U),ti).  (C.9) 

Equation  C.9  will  be  called  the  observation  relation.  Although  the  STK  orbit  data  is 
considered  the  truth  model,  there  is  still  error  within  it.  The  assumptions  made  within 
the  STK  program  create  errors  that  make  the  STK  state  of  the  satellite  differ  from  the 
actual  state  of  the  satellite.  The  difference  between  the  STK  state  and  true  state  can  be 
represented  as  an  error, 

e  =  z  —  z0.  (C.10) 

Applying  Equation  C.9  to  Equation  C.10,  and  then  substituting  in  Equation  C.8  pro¬ 
duces, 

e  =  G(x,t)  —  G(x0,t )  =  G(x0  +  6x,t)  —  G(x0,t).  (C.ll) 

After  expanding  and  subtracting, 

BG 

e  «  ^Sx(t).  (C.12) 

ox 

This  represents  the  true  error  in  the  STK  orbit  data.  It  is  assumed  this  true  error  in  the 
data  will  equal  the  residual  in  the  data,  emr.  The  partial  fraction  matrix  in  Equation 
C.12  relates  the  error  in  the  state  to  the  error  in  the  reference  trajectory, 

Hi(U)  =  (C.13) 

The  residual  can  then  be  calculated  on  the  reference  trajectory, 

r.i  =  zt-  G(xref(ti),ti).  (C.14) 

Remembering  Equation  C.7  and  applying  it  to  Equation  C.14  produces, 

ri  «  Hi5x(ti )  =  t0)5x(t0),  (C.15) 
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r  =  TiSx(t0), 


(C.16) 


where  T  =  HQ.  From  this  relation,  once  the  residuals,  r,  are  determined  the  correction, 
fix  (to),  can  be  found.  However,  the  T  matrix  cannot  simply  be  divided  through  like  a 
scalar  value.  Its  inverse  may  also  not  be  able  to  be  determined,  so  another  method  must 
be  used.  First,  the  transpose  of  T  is  left  multiplied  to  each  side  of  Equation  C.16. 

TTr  =  TTT5x(to).  (C.17) 

The  operation  TtT  is  guaranteed  to  produce  a  square  matrix  which  the  inverse  can  then 
be  taken  to  solve  for  Sx(t0), 

5x  o  =  (TtT)-1Ttt.  (C.18) 

This  change  added  to  the  initial  trajectory  produces  a  better  trajectory, 

X0  =  r  ref  (t0)  +  Sx0.  (C.19) 

This  process  is  continued  until  the  change  is  less  than  a  tolerance  or  a  certain  number 
of  iterations  is  reached. 
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Appendix  D.  Plotting  Results 


D.l  SGP4/STK  Residual  Plots 

In  this  appendix,  all  plots  from  the  discussion  in  Section  4.2  concerning  the  com¬ 
parison  of  the  SGP4  and  STK  propagators  will  be  shown.  As  a  reminder,  coordinate 
residual  plots  show  the  difference  between  the  individual  coordinates  (x,  y,  z)  in  the 
two  propagators.  Residual  magnitude  plots  show  the  magnitude  of  all  the  coordinate 
residuals. 


0  5  10  15  20  25  30 

Time  (days) 


(a) 

Figure  D.l:  SGP4/STK  coordinate  residuals  for  Grace  1  satellite. 
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|  RError  | 
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(a) 

Figure  D.2:  SGP4/STK  residual  magnitude  for  Grace  1  satellite. 


(a) 

Figure  D.3:  SGP4/STK  coordinate  residuals  for  Grace  2  satellite. 
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(a) 

Figure  D.4:  SGP4/STK  residual  magnitude  for  Grace  2  satellite. 


(a) 

Figure  D.5:  SGP4/STK  coordinate  residuals  for  SWIFT  satellite. 
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Time  (days) 


(a) 

Figure  D.6:  SGP4/STK  residual  magnitude  for  SWIFT  satellite. 


(a) 

Figure  D.7:  SGP4/STK  coordinate  residuals  for  Reference  Satellite. 
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(a) 

Figure  D.8:  SGP4/STK  residual  magnitude  for  Reference  Satellite. 
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